{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://cs229.stanford.edu/ps/ps1/ps1.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# (a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use $(k)$ instead of $(i)$ for indexing data point to avoid confusion later, and use $i$ and $j$ to index the dimensions of $x$. Suppose the $x^{(k)}$ is a $n$-dimension vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "J(\\theta) &= \\frac{1}{m} \\sum_{k=1}^{m} \\mathrm{log}\\big(1 + e ^{-y^{(k)} \\theta^T x^{(k)}}\\big) \\\\\n",
    "               &= - \\frac{1}{m} \\sum_{k=1}^{m} \\mathrm{log}\\big(g(z^{(k)})\\big)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $g(z^{(k)}) = \\frac{1}{1 + e^{-z^{(k)}}}$ and $z^{(k)} = y^{(k)} \\theta^T x^{(k)}$. \n",
    "\n",
    "For clarity, in the following derivation, $z = z^{(k)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\frac{\\partial g}{\\partial z} \n",
    "&=\\frac{\\partial \\frac{1}{1 + e^{-z}}}{\\partial z} \\\\\n",
    "&=\\frac{e^{-z}}{(1 + e^{-z})^2} \\\\\n",
    "&=\\frac{1}{1 + e^{-z}} \\cdot \\frac{e^{-z}}{1 + e^{-z}} \\\\\n",
    "&=g(1 - g)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and similarly (unnecessary for this problem, just put here for reference),\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial (1 - g)}{\\partial z} \n",
    "= \\frac{\\partial \\frac{e^{-z}}{1 + e^{-z}}}{\\partial z}\n",
    "= - \\frac{\\partial g}{\\partial z} = g(g - 1)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial J}{\\partial \\theta_i} \n",
    " &= -\\frac{1}{m} \\sum_{k=1}^{m} \\frac{1}{g(z)} \\frac{\\partial g(z)}{\\partial z} \\frac{\\partial z(\\theta_i)}{\\partial \\theta_i}  \\\\\n",
    " &= -\\frac{1}{m} \\sum_{k=1}^{m} \\frac{1}{g(z)} g(z) (1 - g(z)) \\frac{\\partial z}{\\partial \\theta_i}  \\\\\n",
    " &= -\\frac{1}{m} \\sum_{k=1}^{m} \\big(1 - g(z) \\big) \\frac{\\partial z}{\\partial \\theta_i}  \\\\\n",
    " &= -\\frac{1}{m} \\sum_{k=1}^{m} \\big(1 - g(z) \\big) y^{(k)} x_i^{(k)}  \\\\\n",
    " &= \\frac{1}{m} \\sum_{k=1}^{m} \\big(g(z) - 1 \\big) y^{(k)} x_i^{(k)}  \\\\\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then,\n",
    "\n",
    "\\begin{align*}\n",
    "H_{ij} = \\frac{\\partial J}{\\partial \\theta_i \\partial \\theta_j} \n",
    "&= \\frac{1}{m} \\sum_{k=1}^{m} \\frac{\\partial (g(z)  - 1)}{\\partial \\theta_j} y^{(k)} x_i^{(k)} \\\\\n",
    "&= \\frac{1}{m} \\sum_{k=1}^{m} \\frac {\\partial g(z)}{\\partial \\theta_j} y^{(k)} x_i^{(k)} \\\\\n",
    "&= \\frac{1}{m} \\sum_{k=1}^{m} g(z)\\big(1 - g(z)\\big) \\frac{\\partial z}{\\partial \\theta_j} y^{(k)} x_i^{(k)} \\\\\n",
    "&= \\frac{1}{m} \\sum_{k=1}^{m} g(z)\\big(1 - g(z)\\big) y^{(k)} x_j^{(k)} y^{(k)} x_i^{(k)} \\\\\n",
    "&= \\frac{1}{m} \\sum_{k=1}^{m} g(z)\\big(1 - g(z)\\big) x_i^{(k)} x_j^{(k)} \\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly,\n",
    "\n",
    "\\begin{align*}\n",
    "z^THz &= \\sum_{i=1}^{n} \\sum_{j=1}^{n} z_i H_{ij} z_j \\\\\n",
    "           &= \\sum_{i=1}^{n} \\sum_{j=1}^{n} \\big[z_i \\frac{1}{m} \\sum_{k=1}^{m} g(z)\\big(1 - g(z)\\big) x_i^{(k)} x_j^{(k)}\\big] z_j \\\\\n",
    "           &= \\frac{1}{m} \\sum_{k=1}^{m} g(z)\\big(1 - g(z)\\big) \\sum_{i=1}^{n} \\sum_{j=1}^{n} z_i x_i^{(k)} x_j^{(k)} z_j\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that $1 > g(z) > 0$ and $\\sum_{i=1}^{n} \\sum_{j=1}^{n} z_i x_i^{(k)} x_j^{(k)} z_j = ((x^{(k)})^T z)^2 \\ge 0$, \n",
    "\n",
    "$z^T H z \\ge 0$, so $H$ is PSD, i.e. $H \\succeq 0$. \n",
    "\n",
    "Hence, $J$ is convex and it has no local minima."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# (b) Implement Newton's method for logistic regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df_X = pd.read_csv('http://cs229.stanford.edu/ps/ps1/logistic_x.txt', sep='\\ +', header=None, engine='python')\n",
    "# ys = pd.read_csv('http://cs229.stanford.edu/ps/ps1/logistic_y.txt', sep='\\ +', header=None, engine='python')\n",
    "df_X = pd.read_csv('./data/logistic_x.txt', sep='\\ +', header=None, engine='python')\n",
    "ys = pd.read_csv('./data/logistic_y.txt', sep='\\ +', header=None, engine='python')\n",
    "ys = ys.astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_X['label'] = ys[0].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quick look at the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10e470710>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGzdJREFUeJzt3X2MXFd5BvDnib0rTxw2qfAUQhx7UxUFoZQm8SaFRkJq4oQEaNJSVc1KUERWcpCaKrQltCFSEaWgViBCRVGLiV2KiBe1BAQKpUksgigSYK8dJyRxQEB3weFjB1U4ODjEwW//uHfs3c3M7p2Ze+97zrnPTxqtdzy78+7cj/d8H5oZREREzvAOQEREwqCEICIiAJQQREQkp4QgIiIAlBBERCSnhCAiIgCUEEREJKeEICIiAJQQREQkt947gEFs2rTJJicnvcMQEYnKgQMHfmpm7bVeF1VCmJycxNzcnHcYIiJRIblQ5HVqMhIREQABJASS60g+RPJe71hERJrMPSEAuBXAYe8gRESazjUhkNwM4HUA7vKMQ0RE/GsIHwLwDgAnneMQEWk8t4RA8vUAFs3swBqv20FyjuRcp9OpKToRkebxrCFcAeB6kvMAPgXgSpKfXPkiM9tpZlNmNtVurzmMVkREhuSWEMzsdjPbbGaTAG4E8CUze6NXPK46HWD//uyriIgT7z4EmZ0Ftm4Frr46+zo76x2RiDQUzcw7hsKmpqYsqZnKnU6WBI4fP/1cqwUsLABqHhORkpA8YGZTa71ONQRP8/PA+Pjy58bGsudFRGqmhOBpchJ49tnlz504kT0vIlIzJQRP7Tawa1fWTDQxkX3dtUvNRSLiIqrVTpM0PQ1s3541E01OKhmIiBslhBC020oEIuJOTUYiIgJACUFERHJKCCIiAkAJQUREckoIIiICQAlBRERySggiIgJACUFERHJKCCIiAkAJQUREckoIIjI47fKXJLeEQHIDyX0kHyb5GMl3e8UiIgPQLn/J8qwh/BLAlWb22wAuBnAtyVc6xiMia+l0gJmZbJe/o0ezrzMzqikkwi0hWOZY/u1Y/ohnP0+RJtIuf0lz7UMguY7kIQCLAB4ws294xiMia9Auf0lzTQhm9iszuxjAZgCXk7xo5WtI7iA5R3Kuo2qpiC/t8pc0moXRSkPyXQCeNrMP9HvN1NSUzc3N1RiViPTU6WiXv4iQPGBmU2u9znOUUZvkOfm/WwC2A3jCKx4JhIYzxqHdBi67TMkgMZ5NRucCeJDkIwD2I+tDuNcxHvGm4Yyriz1Zxh5/A3iOMnrEzC4xs1eY2UVm9ndesUgANJxxdbEny9jjbwjNVJYwaDhjf7Eny9jjbxAlBAmDhjP2F3uyjD3+BlFCkDBoOGN/sSfL2ONvECUECcf0NLCwAOzdm32dnvaOaG11dJTGnixjj79BgpmHUITmIUhQZmeztvDx8awEvGtXtUks9rH/sccfsaLzEJQQRIbR6WSjZY4fP/1cq5XVbHSzk8AEPzFNJGrqKJUEKSGIDEMdpT40ua1SSggiw1BHaf00ua1y6kMQGYU6SuuhPpuRFO1DWF9HMCLJard1Q6pDt89maULo9tno8y+NmoykmdQWHRf12dRCCUGaR23R8QmtzybRAoX6EKRZ1BYdtxD6bOqekFgCzUMQ6UXzB+LmvTFP4iu3KiFIs6gtWkaReIFCCUGaJbS2aIlL4gUKzz2Vzyf5IMnDJB8jeatXLNIwMa6qKmFIvEDhOQ/hOQB/ZWYHSb4AwAGSD5jZ444xNVcInXV10vwBGdb0NLB9e5LXi+eeyj8ys4P5v38O4DCA87ziaTQNwwxPosMak+HduV2RIPoQSE4CuATAN3r83w6ScyTnOro4ypf4qIkoKUGLE/eEQPIsAPcAeJuZPbXy/81sp5lNmdlUO7FsHITER01ERwlaHLkmBJJjyJLB3Wb2Gc9YGivxURPRUYIWR56jjAhgF4DDZvZBrzgaL/FRE9FRghZHnjWEKwC8CcCVJA/lj9c6xtNcGoYZDiVoceQ27NTMvgqAXu8vK2gYZjgSHtYoYdN+CCIhUoIWB+6jjERkBJqvkLaaj68SgkisNF8hbQ7HV/shiMRI+zqkreTjq/0QJD5q/ihO8xXKF9L553R8lRAkDGr+GIzmK5QrtPPP6fgqIUhxVZWgEl6uobJCp+YrlCfE88/p+CohNMWod6YqS1CJNn9UXujUhMJyhHr+ORxfdSo3waibglfdgVlXB2mNez6ozzciDThY6lSWTBnV4apLUHVUj2tuIw610FmJkDpjh6Hmt1OUEFJXxp2pjg6uKqvHDm3EhT+y2G+moXXGDmu18y/2YzQAJYSuVA96GTfzukpQVe1C5VBcL/SRedxMyzzPQ+yMHUWv8y+VhFeUmUXz2LZtm1Vizx6zVsvs7LOzr3v2VPM+Xrp/38TEaH/f4qLZvn3Z15gsLmZ/N3D60WqZPf545X9P34+sX0xVfrZln+f79mW/a+nfMDGRPZ8Cj2NUEQBzVuAe636TH+RRSUJI6KCvKtabeVlWJsVbbvEtBNR9M63iPE/92kko4RVNCGoyakrvX6Kbghe2tI34wIGs/cazqaPuiUdVnOepd8Y2cPKfEkIDD3pjdZPisWP+hYAiN9My2/urOs9TnguResLrwXtP5d0kF0k+6hZEAw9644VSCFjtZlp2Z2aV53nKtc+UE14PrhPTSL4awDEAnzCzi9Z6faUT02qctCQB6E7WGxvLksGgk/WqVOVEKZ3njVR0Yprrjmlm9hWSk54xnKIdqpol5G0qu+39SxNCt0lr1DiXnOfKDbJS8H0IJHeQnCM514l1fLOEKdSmjhqatJo2vF6KCT4hmNlOM5sys6l2aBeuCDBc5+9qP1Nxv1Zq88mkPMEnBJGgDVPULvIzFXZmNmWktQzOfbXTvA/hXvdO5YLU7iqnDNP5G8DKmgGEIDWLYrVTkrMAvgbgQpJHSM54xrMWtbvKMsMUtQMonmuktfTjXkMYhGcNQaUqeZ5IawhLQ1FttxmiqCHEJICCnYRmmKJ2QMXzUAdZiR/VEAoKqGAnoRmmqK3iudQoiolpMekW7FZObq36WtZ9IwLDTGrUREgJkJqMBlD3sibqxBaROqnJKFDdJqqNxzuYxDzmMYmnW201UYnIwNSpHLn5eWAas1jAVjyAq7GArfgTm1UntiQn1d1rY6SEEKgLzurgw8dncCaO4xwcxZk4jo88M4MLzmrIVaO7RCOoWTQsSgiB2nRsHutby8e5rm+NYdOxeZ+A6qS7RDxGSNxaUyk8SgihmpzEOJaveDmOBuzkprtEPEZM3JrbEx4lhFAFNIGpVrpLhGtpbaCExB3KxnVymhJCyBq2fR8A3SVCtbI28NGPjpy4Ry7zlN3PpH4rJYTgNW19gabWjELWqzbw3veWkriHLvOU3c+kfisAmoeQrE4HeOih7N+XXBLh/VRTtMOxf392ozx69PRzExPAbbcB73tf/ftSl72OTAPWpdHSFQ02Owu8+c3ZNQpkNfuPfzyyFict7RCOfs14N9+cPepO3GXvOV3lHtaRUZNRYjod4KabTicDILuWe/X3qck0Qh4HbbVmPI8mzbL7mdRvdYoSQmLm54F1657//BlnLO/vU5NphJYcNNu6Fd/9+9n68kJIAxzK7mdSv9Uprn0IJK8F8E8A1gG4y8z+YbXXqw9hbZ0OsGUL8Mwzy59f2iQaQpOpuggG1OOg/QItvGzDAv5xdzuu5sCylH0SJXxSBr+WEcl1AD4C4DoALwcwTfLlVbxXk5pG2m1g9+6sCbRrfHx5gcd7qL9qJ0PocdBOYAwvemY+rXl7g1ysZTdXNW1EXw+eTUaXA/iOmX3PzJ4F8CkAN5T9Jk28+UxPA08+Cdx3X/Y4cmR5Dd+zyVQTkYfU46CN4QTmMZnOvL0mXqyB8UwI5wH4wZLvj+TPLUNyB8k5knOdAe8aTb75tNvANddkj5UFHs8mU+/aSVd0tcb8oFmrhaOYwC/Qwk3YhZ+inUb/Z10Xa3QHvl6eCYE9nnteh4aZ7TSzKTObag94xwrl5hMirz7CEAZ0RFsQnZ4GFxbw9ffsxcs2LOCLE9Pp9H+OerEWudFHe+BrZGZDPQC8ZdifzX/+VQDuW/L97QBuX+1ntm3bZoNYXDRrtcyA049WK3te/OzZkx2HiYns65499b13KufE4qLZvn3xxd3XKAeme0KdfXb/EyqVAz8kAHNW4L48Sg3h3SPmov0AXkryApLjAG4E8PkRf+cyGk0WJs8RjKnUGqPp/yxScu+O7rnzzsEv1qJNTakc+IqtOlOZ5CP9/gvAi0Z5YzN7juQtAO5DNux0t5k9Nsrv7GV6Gti+PdnRZNHymogcQpNVY8zOZjfn8fHsQ++1tMXK19x5J3DppcUv1qKzjEc98AkPSV1mteoDgJ8AuBjA1hWPSQA/LFIFKfMxaJORSC+eTVaNUaSJpoxmnEF+x7AHvkiTVOBQsMlorbWM7gVwlpkdWvkfJL9cXloSqY9qjTUoUnIvYw2hbrvwzMzyRfZ6/fwwB35pk1Q3zpkZ4IUvjHTVyNVptVMRKV+R6fBlTpmvqkmn10qvALBxI3DyZH0rvI4o+JnKIpKwIiM6yhz1UVUve6++BwB4+ukkJzaphiAi1SlScg+9w7bb8X3GGVkiWGpiIhsud9llPrEVpP0QRMRfkeFkoe990e17eOgh4IYblq8cmdgQNTUZ9aEZ7iJySnctmN27k57YpITQg2a4i0hPIe0LUQH1IawQwl4BRRRtdg29eVZEqqdRRkOKYYZ70RqMajoiMgjVEFYIvYZQNL7Q/w5pCFVRg6AawpDqWhBv2E7rojWYGGo6kjhVUaOjhNBD1f1Go1wnRdfo0iJuvWn0WE2avDtVxJQQ+qhq4uOo10nRGoyW/n4+FVhrFFsVVSUFAEoItZ8HZVwnRWswpdd0Ir5oVGCtWUxVVJUUTml0QvA4D8q6TorWYEqr6QzwYYWYN2IrsEYvliqqSgrLNDYheJ0HsVwnywzwYYVa2IqpwJqMGCZxqaSwjEtCIPnHJB8jeZLkmkOhquB5HsRwnSxT8MMKubAVZSJOQeh7faqksIxXDeFRAG8A8BWn93c/D0K/TpYp+GGFXtjySMQhNp/JEsOWFBI9sC4JwcwOm9m3PN67SyXGART8sLyTbBF1JuJQm89khUFLCgkfWNeZyvk2nG83s77Tj0nuALADALZs2bJtYWGh1Bg0kXIABT6s7tLxS3czDL5JrAKaKZ6oSA+s+34IJPcCeHGP/7rDzD5X9PeY2U4AO4Fs6YqSwjsl9KXYg1Lgw9J+xZkytguWACV+YCtLCGa2varfLWFTko2j+UyGkPiBbeywU5EqqY8qUYkfWJc+BJJ/CODDANoAfgbgkJm9Zq2f057KEhv1USUqsgPr3oewGjP7LIDPery3SJ3UfJaoRA+smowGMMrQ40SHLYtIQpQQChpl6HHCw5ZFJCHaMa2AUYYeRzpsWUQSoh3TSjTKkgyhL+cgItKlhFDA5CTwy18uf67o0OPEhy1LYNRXJaNQQihg717g5MnT34+NFR96nPiw5SA19aaovioZlfoQ1lBWH0Bkw5aj1V1LaXw8q5k1ZS0l9VXJatSHUJKy+gDqWmWzqtJxnaXuYd8r5P0YqhZKX1VTa2epUEJYQ0x9AFU1GdTZFDHKe4VyU/QQwnmqJqsEmFk0j23btpmHPXvMWi2ziYns6549LmGsanExiw04/Wi1sudD/L1VvFedsYbI8zxt+mcfOgBzVuAeqxpCATFseVlV6bjOUveo79X0DnzP87TJtbOhBdi+5rKWUYxCX7qkqiaDOpsiynivpu/H4HWehtBkFZVARz+ohpCIqkrHdZa6h3mvXoWsqParTkTTa2cDCXj0g4adJqaq4a1l/N6iv6Po6wItZDWahlcXsH9/1vN+9Ojp5yYmsra+yy6r5C2LDjtVQhiATvbhlX3zXnXcPXSgJGAOk0Y0D6FkGlI3vCpqyP06MZ/6qA6UBC7g9jWvHdPeD+D3ATwL4LsA3mJmP1vr52Jc7VSqqSH3Oibnb+hggVtBHSiJQY1NDqHXEB4AcJGZvQLAtwHc7hRHIRpSN5oqRqD0KmR97I55UAdKYhHg6AeXhGBm95vZc/m3Xwew2SOOomIaUhfg0ObKasgrx92/5ubJeA6USIBC6EO4CcAX+/0nyR0k50jOdZzucgE3+S0Tcj9HVZOmlhWyKjhQISbYaOjDi05lfQgk9wJ4cY//usPMPpe/5g4AUwDeYAUC0Sij/or0c4Qcf6lK+kM1rHUE+vCCEvywU5JvBvBWAFeZ2S+K/Ix3QgjZWh23uj4Ho4EEI9CHF5ygO5VJXgvgrwFcXzQZyOpW6+cIeGJksDSQYAT68KLl1YfwzwBeAOABkodI/qtTHMlYrflc1+fgYhpIEBx9eNFyWdzOzH7T431T129hN12fg+sm2JmZLHmeOBHmQIIg6cOLlpauaIhuH8LS61N9CGtrTEd8FfThBSP4TuVhKCGMRtenL33+4iXoTuUYpDiEOsCJkY0R8hwRkS4lhB5Su3hTTG4x0SgviYUSwgqpXbypJbcYaZSXxEIJYYWULt7UklusNMpLYqGEsEJKF29KyS1msayFJaKEsEJKF29KyS12VS3uJ1Iml4lpoes3wSs2mh8UlnZbn72ETQmhj1QuXq/kpjH38dMxbB41GTVA3fMPNLIpfjqGzaSZypEKtfSmlY/jp2OYHs1UTljIpTfvkU1ek/BSmvznfQzFjxJCZEKfW+A5sskrUYacoIeh0WnNpYQQmdBLb17Ddr0SZegJehixD71OqbZWNyWEyMRQevMYc++VKENP0MOKdd5EarW1unltofkeko/ku6XdT/IlHnHEKJbSW90jm7wSZQwJelixrY6bYm2tbl41hPeb2SvM7GIA9wL4W6c4ohRq6c2zql57osz/2DY6USToJki1tlYnry00n1ry7UYA8Yx9DURoE+e6O7KNj2clZo8d2WqbhLfij53etQvbF6aDHAbcJCnX1uriNg+B5HsB/CmAowB+z8x6litJ7gCwAwC2bNmybWFhob4gpZBGjVtv1B8bH20V25v7PASSe0k+2uNxAwCY2R1mdj6AuwHc0u/3mNlOM5sys6m2LrggNaqq3qg/Nj6hNqfGorImIzPbXvClewB8AcC7qopFqtWoqnqj/tg4hdacGhOvUUYvXfLt9QCe8IhDyhHLyKdSNOqPlaZx6UMgeQ+ACwGcBLAA4K1m9uRaP6e1jMIW6vpKlWjUHyuxK9qH4DXK6I883leq1aiqeqP+WGkKzVQWEREASggi0dFaPVIVJQSRiGitHqmSEoIEKbZScB3xaq0eqZoSggQntlJwXfFqTpxUTVtoSlBiWxmiznhj+2wkHO5LV4gMI7ZScJ3xak6cVM1lHoJIP7GtDFF3vLWt6CqNpBqCBOed7wQ2bIijFOxRao9t4xqJh2oIEoyl2wyQwG23ATffHP6NT6V2SYU6lSUI6jAVqY46lSUqsXUmi6RICUGCEFtnskiKlBAkCBpSKeJPncoSDHXOivhSQpCgaJsBET+uTUYk307SSG7yjENERBwTAsnzAVwN4PteMUj9YlvFVKRJPGsIdwJ4B4B4JkLISGJbxVSkaVwSAsnrATxpZg8XeO0OknMk5zoqVkZLa/mLhK+yTmWSewG8uMd/3QHgnQCuKfJ7zGwngJ1ANlO5tAClVt2JZ0tnIncnnqkTWSQMlSUEM9ve63mSvwXgAgAPkwSAzQAOkrzczH5cVTziSxPPRMJXe5ORmX3TzH7dzCbNbBLAEQCXKhmkTRPPRMKneQhSG008Ewmbe0LIawnSEJp4JhIurWUkIiIAlBBERCSnhCAiIgCUEEREJKeEICIiACLbU5lkB8CCdxwANgH4qXcQBSjO8sUSq+IsXyyx9opzq5mtOb4vqoQQCpJzRTas9qY4yxdLrIqzfLHEOkqcajISEREASggiIpJTQhjOTu8AClKc5YslVsVZvlhiHTpO9SGIiAgA1RBERCSnhDAAkrtJLpJ81DuW1ZA8n+SDJA+TfIzkrd4x9UJyA8l9JB/O43y3d0yrIbmO5EMk7/WOZTUk50l+k+QhknPe8fRD8hySnyb5RH6uvso7ppVIXph/jt3HUyTf5h1XPyT/Ir+WHiU5S3LDQD+vJqPiSL4awDEAnzCzi7zj6YfkuQDONbODJF8A4ACAPzCzx51DW4bZDkkbzewYyTEAXwVwq5l93Tm0nkj+JYApABNm9nrvePohOQ9gysyCHjNP8t8B/I+Z3UVyHMCZZvYz77j6IbkOwJMAfsfMQpgPtQzJ85BdQy83s+Mk/wPAf5nZx4v+DtUQBmBmXwHwf95xrMXMfmRmB/N//xzAYQDn+Ub1fJY5ln87lj+CLKGQ3AzgdQDu8o4lBSQnALwawC4AMLNnQ04GuasAfDfEZLDEegAtkusBnAngh4P8sBJC4khOArgEwDd8I+ktb4Y5BGARwANmFmScAD4E4B0ATnoHUoABuJ/kAZI7vIPp4zcAdAD8W94MdxfJjd5BreFGALPeQfRjZk8C+ACA7wP4EYCjZnb/IL9DCSFhJM8CcA+At5nZU97x9GJmvzKzi5HtrX05yeCa4ki+HsCimR3wjqWgK8zsUgDXAfizvKkzNOsBXArgX8zsEgBPA/gb35D6y5u0rgfwn96x9EPy1wDcgGzP+pcA2EjyjYP8DiWEROVt8vcAuNvMPuMdz1ry5oIvA7jWOZRergBwfd42/ykAV5L8pG9I/ZnZD/OviwA+C+By34h6OgLgyJIa4aeRJYhQXQfgoJn9xDuQVWwH8L9m1jGzEwA+A+B3B/kFSggJyjtrdwE4bGYf9I6nH5Jtkufk/24hO6Gf8I3q+czsdjPbnG/3eiOAL5nZQCWvupDcmA8kQN4Ecw2A4EbFmdmPAfyA5IX5U1cBCGrQwwrTCLi5KPd9AK8keWZ+D7gKWf9hYUoIAyA5C+BrAC4keYTkjHdMfVwB4E3ISrLd4XKv9Q6qh3MBPEjyEQD7kfUhBD2kMwIvAvBVkg8D2AfgC2b2384x9fPnAO7Oj//FAN7nHE9PJM8EcDWyEnew8trWpwEcBPBNZPf3gWYta9ipiIgAUA1BRERySggiIgJACUFERHJKCCIiAkAJQUREckoIIiMieS3Jb5H8DslgZ9uKrEXDTkVGkK+A+W1k49SPIJtPMR3ayrIiRaiGIDKaywF8x8y+Z2bPIlva4gbnmESGooQgMprzAPxgyfdHEOBS4yJFKCGIjIY9nlM7rERJCUFkNEcAnL/k+80YcFMSkVAoIYiMZj+Al5K8IF8z/0YAn3eOSWQo670DEImZmT1H8hYA9wFYB2C3mT3mHJbIUDTsVEREAKjJSEREckoIIiICQAlBRERySggiIgJACUFERHJKCCIiAkAJQUREckoIIiICAPh/H1nWqC865OUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.axes()\n",
    "\n",
    "df_X.query('label == -1').plot.scatter(x=0, y=1, ax=ax, color='blue')\n",
    "df_X.query('label == 1').plot.scatter(x=0, y=1, ax=ax, color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Newton's method addresses getting to $f(\\theta) = 0$, and minimizing $J(\\theta)$ means getting $\\frac{\\partial J}{\\partial \\theta}$ to 0. There after applying Newton's method, extending it to multidimensional setting (Newton-Raphson method), the update rule becomes:\n",
    "\n",
    "\\begin{align*}\n",
    "\\theta &:= \\theta - \\frac{\\partial J(\\theta) / \\partial \\theta} {H} \\\\\n",
    "       &:= \\theta - \\frac{\\nabla_{\\theta} J(\\theta)} {H} \\\\\n",
    "       &:= \\theta - H^{-1} \\nabla_{\\theta} J(\\theta)\n",
    "\\end{align*}\n",
    "\n",
    "Note, the $H$ in the denominator may not be a valid mathematical expression as it is actually an inverse operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xs = df_X[[0, 1]].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(99, 2)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xs.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adding a columns of ones for the intercept terms, and also use column vectors\n",
    "Xs = np.hstack([np.ones((Xs.shape[0], 1)), Xs])\n",
    "ys = df_X['label'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copied equations from above:\n",
    "\n",
    "$$z^{(k)} = y^{(k)} \\theta^T x^{(k)}$$\n",
    "\n",
    "$$g(z^{(k)}) = \\frac{1}{1 + e^{-z^{(k)}}}$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial J}{\\partial \\theta_i} \n",
    " = \\frac{1}{m} \\sum_{k=1}^{m} \\big(g(z^{(k)}) - 1\\big) y^{(k)} x_i^{(k)}\n",
    "$$\n",
    "\n",
    "$$\n",
    "H_{ij} = \\frac{1}{m} \\sum_{k=1}^{m} g(z^{(k)})\\big(1 - g(z^{(k)})\\big) x_i^{(k)} x_j^{(k)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "converged after 7 iterations\n"
     ]
    }
   ],
   "source": [
    "all_thetas = [] # collect for demonstration purpose\n",
    "theta = np.zeros(Xs.shape[1])\n",
    "tol = 1e9\n",
    "n_iters = 0\n",
    "while tol > 1e-6:\n",
    "    zs = ys * Xs.dot(theta)\n",
    "    gzs = 1 / (1 + np.exp(-zs))\n",
    "    nabla = np.mean((gzs - 1) * ys * Xs.T, axis=1)\n",
    "    \n",
    "    # Refactor, more efficient way of calculating hessian\n",
    "    hessian = np.zeros((Xs.shape[1], Xs.shape[1]))\n",
    "    for i in range(hessian.shape[0]):\n",
    "        for j in range(hessian.shape[0]):\n",
    "            if i <= j:\n",
    "                hessian[i][j] = np.mean(gzs * (1 - gzs) * Xs[:,i] * Xs[:,j])\n",
    "                if i != j:\n",
    "                    hessian[j][i] = hessian[i][j]\n",
    "            \n",
    "    delta = np.linalg.inv(hessian).dot(nabla)\n",
    "    old_theta = theta.copy()\n",
    "    theta -= delta\n",
    "    all_thetas.append(theta.copy())\n",
    "    n_iters += 1\n",
    "    tol = np.sum(np.abs(theta - old_theta))\n",
    "print('converged after {0} iterations'.format(n_iters))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10abfaf28>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAEKCAYAAABNOm93AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW9+P/XZzLZyZ4QQgJJ2BLCGpYgRVzY3FhEa1tsba/FqtelWlt7a7231W9/dnNt7eL1itW2lrqiiLiAKyAgSwABE9YECEsSCNm3yXx+f0zSBswymczMOWfm/Xw8eIRMZua8Z87MeX/2j9JaI4QQQoi+sxkdgBBCCGFVkkSFEEIID0kSFUIIITwkSVQIIYTwkCRRIYQQwkOSRIUQQggPSRIVQgghPCRJVAghhPCQJFEhhBDCQ3ajA+iL5ORknZWVZXQYQghhKdu2bavUWqcYHUcgslQSzcrKYuvWrUaHIYQQlqKUKjU6hkAlzblCCCGEhwxPokqpEKVUoVJqldGxCCGEEH1heBIF7gK+MDoIIYQQoq8M7RNVSmUAVwEPAfcYGYsQQoiebdu2baDdbn8GGIs5KmH+4AR2OxyOmyZPnlx+/h+NHlj0BPBjIMbgOIQQQvTCbrc/M2jQoNEpKSlVNpstKDajdjqdqqKiIu/kyZPPAAvP/7thJQml1HygXGu9rZf73ayU2qqU2lpRUeGn6IQQQnRhbEpKSk2wJFAAm82mU1JSqnHVvr/8dz/H09kMYKFSqgT4JzBLKfX38++ktX5aaz1Faz0lJUWmOQkhhIFswZRAO7S/5i7zpWFJVGt9n9Y6Q2udBXwD+EBr/S2j4jFURQVs2eL6KYQQwjKCpWPYvJYvh8xMmDvX9XP5cqMjEkII08rPz88FKC4uDnvqqacS+/t8d955Z/qgQYPGR0VF5XvyeFMkUa31R1rr+UbH4XcVFbB0KTQ2QnW16+fSpVIjFUKIbhQWFhYB7N+/P/zFF1/sUxJ1OBxfuu3qq68+u3nzZo+nWZoiiQatkhIICzv3ttBQ1+1CCCG+pKPGeP/996dv3bp1QG5ubt6DDz440OFwcMstt2SMHTt29KhRo/IefvjhZIBVq1bFTJs2bdSCBQuyc3Jyxpz/fLNnz67PzMxs9TQeo6e4BLesLGhpOfe21lbX7UIIIbr10EMPlT366KOpH3744QGARx55JDkuLq5t9+7dXzQ2NqqpU6fmLliwoAZg165d0YWFhXtyc3Nben7WvpMkaqSUFFi2zNWEGxrqSqDLlrluF0IIk/vLhsNJx6oaw731fBkJkc03zsg+7clj165dG1tUVBS1cuXKBIDa2tqQvXv3RoSFhenx48fX+yKBgiRR4y1ZAnPmuJpws7IkgQohLMPThOcLWmv16KOPHrn22mtrOt++atWqmKioKKevjit9omaQkgJTp0oCFUIIN8XFxbXV1dWFdPw+d+7c6j//+c8pzc3NCmDXrl3hNTU1Ps9xkkSFEEJYTkFBQaPdbtc5OTl5Dz744MAf/OAHlbm5uU3jxo0bPXLkyDHf+973MltbW1Vvz3PrrbdmpKamjm9qarKlpqaOv+eeewb3JQ6ltXUWn5gyZYqWTbmFEKJvlFLbtNZT+vs8O3fuLJkwYUKlN2Kymp07dyZPmDAh6/zbpSYqhBBCeEiSqBBCCOEhSaJCCCGEhySJCiGEEB6SJCqEEEJ4SJKoEEII4SFJokIIISzDm1uh1dbW2i655JIR2dnZY0aMGDHmtttuS+/rc0gSFUIIYRne3grthz/84anDhw/v2b17997NmzcPeOmll2L78pySRIUQfVdRAVu2yN63wu+8uRVaTEyMc8GCBbUAERERevz48Q1Hjx4N+/JRu2fYAvRKqQjgEyC8PY5XtNY/NyoeIYSbli937TwUFubaym/ZMtdGCkL4kbe3QqusrAxZs2ZN/L333nuqL3EYuYtLMzBLa12nlAoF1iul3tZabzIwJiFETyoqXAm0sdH1D1y/z5kjGygEo81PJVFV6rWt0EjIbGbarX7fCq21tZVrrrlm2M0333wqLy+vT1umGZZEtWvR3rr2X0Pb/1lnIV8hglFJiasG2pFAwbUXbkmJJNFg5GHC84X+bIV2/fXXZw0bNqzpZz/7WXlfj2ton6hSKkQptQMoB9ZorTcbGY8QohdZWa4m3M5aW123C+FH3toK7fvf//7gmpqakGXLlh31JA5Dk6jWuk1rPRHIAAqUUmPPv49S6mal1Fal1NYKGcQghLFSUlx9oJGREBvr+rlsmdRChd95Yyu0gwcPhj755JNp+/fvjxgzZkxebm5u3mOPPZbclzhMsxWaUurnQL3W+pHu7iNboQlhEhUVribcrCxJoBYgW6H1n+m2QlNKpSil4tv/HwnMAYqMikeYhEydsIaUFJg6VRKoCHpGNuemAR8qpXYBW3D1ia4yMB5htOXLITMT5s51/Vy+3OiIzMXqBQyrxy9EFwxLolrrXVrrfK31eK31WK31/zMqFmECnadOVFe7fi5dKhfcDlYvYFg9fiG6ISsWCXPomDrRWcfUiWBn9QKG1eMXogeSRIU5yNSJ7lm9gGH1+IXogSRRYQ4ydaJ7Vi9gWD1+IXogSVSYx5IlUFoKa9e6flphPVZ/DJaxegHD6vELU/HmVmgAM2fOHJmTk5M3YsSIMddff/3QrnZ66YkkUWEuVpo64c/BMlYsYHRm9fiFaXh7K7Q33njjYHFx8d59+/btOX36dOizzz6b0JfnlCQqhCeMGCxjpQJGV6wevzAFb26FBpCYmOgEaG1tVa2trUqpHhc5+hIjd3ERwrpkIXYhDOXNrdAuvPDCkbt27Yq++OKLq2+88caqvsQhSVQIT8hgGWPIcoOm8sIXLySV1ZV5bSu09AHpzd8c/U2/b4W2fv36/Q0NDWrx4sXD3nzzzdjFixfXdHff80kSFcITHYNlli511UBbW2WwjK/JZuCm42nC84X+bIUGEBUVpefPn392xYoV8X1JotInKoSnZLCM/8iCDeI83tgKrbq62lZaWhoKro2533nnnbjc3NzGnh5zPqmJCtEfKSlS+/QH6YMW5+m8Fdr1119f+d///d/lJSUl4ePGjRuttVaJiYmtq1evPtjTc9TU1NiuuuqqES0tLcrpdKoZM2bU3HvvvX0qmZlmKzR3yFZowmukb81aKipc04g6J9HISFcLgJy/XslWaP1nuq3QhDCMLIZuPWZbsEF2pBHtJImK4CJ9a9Zllj5oKYSJTiSJiuAii6Fbm9ELNkghTJxHkqgILjK/U/SHFMLEeSSJiuBitr41YS1SCBPnMSyJKqWGKKU+VEp9oZTao5S6y6hYRJAxS9+asB4phInzGFkTdQA/1FqPBi4AbldK5RkYT3ALttGGRvetCeuSQpihvL0VWodZs2aNGDly5JcWqO+NYUlUa31Ca729/f+1wBdAuq+O99yGw7y89ShNrW2+OoR1yWhD8wm2Qo3VSCHMMN7eCg3g+eefj4+OjvYoOZiiT1QplQXkA5u7+NvNSqmtSqmtFf24oPzHjGzGZ8Tzpw8P8Lu1+ymprPf4uQKKjDY0HynUCNEtb2+FVl1dbfv973+f+sADD5zwJB7Dl/1TSg0AXgXu1lp/adFfrfXTwNPgWrHI4wOd2kPOwNHkzMuhrtnBGzvK+MdnR5ialcis3IGE2Pq2h1zAkOXUzKVzoabjnCxdCnPmyPkQohNvbYV2zz33pN91112nBgwY0Osi9V0xNIkqpUJxJdAXtNav+fJYzvpKbLsehMgEBoz/Gt+clonWmq2lVTzyXjEJUaEszs8gJcZru/pYg4w2NBcp1AgLOfO3vyW1HvPeVmihGenNiTfc4Let0D799NPIw4cPhy9btuxocXFx2JeftXeGJdH27cOXAV9orR/z9fHef3I39eGxDLlhDJN2r0DVnURlzmDqyLlMzUqkoraZFYXHOFPfyuzRA5mSmUBfdzi3JNnSy1ykUCMsxNOE5wuebIW2bt26Abt3745KT08f53A41JkzZ+wFBQU5n332WbG7xzWyT3QGcAMwSym1o/3flb462KUP30LmkGHU/W8JHyxXvHowhTNKwwe/gHWPkdJWzs0XDedH80ZxtqGVX79dxAubS6lr7rojOqDIaEPzkCkUQrjFG1uh/dd//VdFeXn5rrKyss8/+eSToqysrOa+JFAwsCaqtV4P+K2qZ7fbyb99MQDVJSfZ+dRqPn/yKPX2eCIWp3Fx0SpCao5jHzKNuTmXMTcvlZLKepatO0yb08n8CYMZlRrjr3D9T7b0Mo8lS1x9oLLLjBDd8sZWaN4Q9Fuh7V+xjuMffYHTqakaWMuka0aSdWIHhEXBuOsgIYum1jbe3HmcA+V1jE2P47Ixgwizm2JgsxBC9Eq2Quu/7rZCM3x0rtFGLp7JyMUzaaqpp/DJFZQ+Vc4+lUjLjDDmhL1NRE0ZEelTuG7ilWAfws6jZ3li7T6iw+0szk9ncHyk0S9BBDPZFzWwyfk1vaBPoh0iYqOZfv+3ADix5Qv2LV/HpvUh1EQnkrW4gfEf/xrsEUwYey0TLs/lbEMLr20v41RNExeOTGbG8GRswTpNRhhj+XLXgLCwMNdgpGXLpD87kMj5tYSgb87ticPh4PNn3qJm7ykcuo2aUU3MmhpOXM1xSJsIoxfgtIWx/kAlGw5UMigugmvyM4iLCvVbjCJIVVS4FmLoPBUmMtI1MExqLNbn5fMrzbn9J825HrDb7eTfugiA2mPlFP7xTXa80EZ9SAoxl51hRvlvsYWEctGYa7joytGUnW3k75tLaWhxcPmYNMZlxBn8CixGmq7cJ/NJvc9Mnz85v5YhSdRNMRkDuehXSwHYv3I9x9/Zy8fOgZxNqmNqyDtkFP6V9NRx3H7hQlpUOO/uOcmqXccZmRrD/PFpRISG9HKEICdNV30j80m9y2yfPzm/liHNuf3QUt/E1t+9SuuJeppVC21THMxOryHMZoe8qyE1j+KTtby16zghNhtX5w8mMyna6LA956uSegA3Tfq0ctNx4e+8SIYUPPrOrJ8/L55fac7tP2nO9YGw6Ai+8tNvAnByezFFL3zMp1uiqI2oJbt1NWP5Bzkpo8m5ZDF1OozXC8t4YfMRCrISudTf6/X292ruy5J6gDZd+bxyI/NJvcOsnz85v13Kz8/PLSwsLCouLg778MMPB9x6661n+vN8BQUFOeXl5aERERFOgPfff39fenq626vsSE3UyxwOB7ufe5uzu07Qptuoy6pl1sRGYlQIjF6IHjSOLSVVfFBUTlJ0GIsnpZM8wMfr9fb3au7rkrq/agJ+7PMya+VGdCEITlYg1kRXrVoV03kBenc4HA7s9nPrjgUFBTmPPPLI0Ysuuqihp8d2VxOVFQO8zG63M/GmBVzy+5spuO86Ek4nsf31NN5eGcOn696C9/6HgjNv8pNZ6SzKH8yr247xm3eK2FpyBp8UaLyx1VlHSb2zjpK6N/hjqTs/by/m67fMVKy+96kstWgp3t4Krb+kOdeHYgYnc9EvXYORDr69kSPv7uKjtkRq4k8ytfy3DA4L5Zbc+ThSJ/BBcQW/fruIzKRoFk4czIBwL50abzRV+WOQgy+brgzYXsztt8xMI0I9YbYBOZ7q6fNn9XMUoLy1FRrATTfdlGWz2ViwYEHVb37zmxM2m/v1S0miHXz8RRl+xXSGXzGdlqZmtv3uZfavbmAPLTjHvMms7FeZlziMebO+yuE6O8+sO4TTqb2zXq83EqC/dnrx1fq9BvR5ufWWGZGAvPk5D7S9T7v6/AVKIcFHdn1wNKnmdJPX+qNikyKax88a4ret0ABefPHFQ9nZ2a1VVVW2+fPnD//Tn/6UdMcdd7gdgyRR8OsXJSwinOn/5VoZqWLnAXb/9X027FHUhVUx/MhvGR0Xxt05V9I4MJ83Pz/BK9uOMT4jjnl5Hq7X660EaOVBDt0VJAYMcDVD+uj19PiWGZGAvP05N+uAHG8JtEKCD3ia8HzBk63QALKzs1sBEhISnF//+tfPfPbZZ9GAJFG3GfhFSZkwgksfHUFbWxt7/vo2pz4u44R2Uj9oNbMmvcrXkobDpV9lZ6X613q910xKJy2uj+v1eisBWnWnl64KEkuXwuTJPi84dfuW+TsB+eJzHuhzGQO9kGBx3W2FNn/+/Nrw8HC9a9eu8KysrNaenqO1tZXKykp7Wlqao7m5Wa1evTpu1qxZtX2JQ5KoCb4oISEhjL9xPtwI9RVVbH3iNbaubqPJVk/c+F8zbWgEE0ZfTlXSJFbsOM6p2iZmjkjhK8OT3F+v16oJ0Fs6FyQGDHAlUCNrGP5OQL74nAf6hu6BXkiwOG9shdbY2GibM2fOyNbWVuV0OtXMmTNr7rnnnj6NkJMpLiYe3n7wvc0cWb0D2qB+QCWTL6giLXUEzrHXsf64kw0HK0mLjWCxrNfbN1u2uEbpVlf/+7bYWNem5FOn+i+O3ibTe7v/0lef80AeeBMgC1oE4hQXfzPlYgtKqWeB+UC51nqsIUGYuDQ9fN40hs+bhqOphS1Pvkzx2hT20Ioe9jCzxoRxUe4cyuKn/Gu93ivGpjE2Xdbr7ZVZahg9NbN7u//Sl5/zQG7lsPJYAOEXhtZElVIXAXXAX91Joj5dbMEipenTew+zc9l72FoUDfYqRhWcZETWCFryvso7JW3sKatmVGoMV8l6vT0zcw1Dao3Cy6Qm2n+mrIlqrT9RSmUZGcO/WKQ0nZSXzaxHbwFg5/NvcXRjPEc2OGlM/h2zCjQLc+ZSFJXGnz48gD3ExqKJFl+v11fMXMPwZT99p8+55FMh+s/0A4uUUjcDNwMMHTrU4GjMZcJ3roLvQENlNZ/97mW2ve2k6e1dJI5byQ/ysqnPvYbX91fyj81HKMhO5JIcP6/Xa3ZmLTj5oblZpj8K4R2GDyxqr4muMrw5N0Ac/mALh9/Yhq1NURd5kmkz60jOnctntol8uP+M/9brDSaeVOl6e4wPm5tNPJZO+Ig05/afKZtzhfdlz5pK9qypOJpa+OyPL7NnTR2O94qxDX2XH09P58zwq3ll2zHONrQyZ/RAJmcmoJTUTj3mSZXOncf4sLnZBLO6hAgYsgB9H1llrW17RBhf+eE3ueTJW8i/cxGUj+STF6PZ8puXmHfsT9ybdZAztQ386u0i/rH5CPXNbu/8Izp4srh/Xx6TkuKacuPlzGaWwclCeCI/Pz8XoLi4OOypp55K7O/zNTU1qSVLlmRmZWWNzc7OHvPcc8/F9+XxhiZRpdRyYCOQo5Q6ppRaamQ8vfHzRiBek5QzlFmP3Molf7iFwQUTKN08ik/+eATni3/g7og3uDC5gac/OcRj7xWz/1SfFusIbp5s1WKC7V1k0xJhZYWFhUUA+/fvD3/xxRf7lEQdji9XFu677760lJSU1pKSkt0HDhzYc9lll9X15TkN7xPtCyP7RAOtH6nxzFk2P/EKVLXRoutJHnOMMRdezMrGceyvbGZCRjzzxqQSGiKNFd3y5ENhog+SjM4NHoHUJxoVFZXf0NBQOGHChNxDhw5FpKentyxZsqTy/vvvL7/99tszNmzYENPS0qK+973vld97772Vq1ativnFL36RNnDgwNa9e/dGHTx4cE/n5xs0aND4ffv27Y6Nje12fV2QPtF+C7R+pMjEeC75fzcBcPjjrRxesYWNe04QFV7ILXNCKA9ZxBNrq4kK83C93mDgyQIGJlrcw6yDk4Vwhze2QqusrAwBuOeeewZ/+umnMZmZmc1PP/30kSFDhrjdvyVJ1E2B3I+UffEUsi+egsPhYMuTL7FrdS1t+mOmZpQyeW4Br+50cqK+jYtGutbrlYFInXgyAMjMc1SF6IPtb69Mqik/5b2t0AamNk+6YqHftkJrbW1Vp06dCr3wwgvrnnnmmWMPPPBA6p133jnk9ddfP+zucSWJusmoCoQ/m9zsdjvTf3A9AGcOHKPwqTfZ9pczpNpe46LpdVSeXciv36lgcFwkV+enExcp6/UCnlXppBooAoCnCc8XPNkKLTU11REREeG84YYbzgJ861vfOvP3v/89uS/HlQ6vPliyxNV1tXat66evJ6cbOZApcUQGsx/5Ty75w60Mnj6B0o3p1D+3lcnr/sacto/4x4ZiHn63iN1l1b0/mRBCeFl3W6E1NzcrgF27doXX1NT0mONsNhuzZ8+ufuutt2IAVq9eHTty5MjGnh5zPqmJ9pG/KhAdMyGiGyvIaiyhhCyWLk0xZD/gCUuugCXQeLaWTY//k/0vNjGC90jKq+BU1BW8uTOFnEExXDlO1usVQviHN7ZCA3jssceOXX/99dk/+tGPQpKSkhx//etfS/oSh4zONaktW+Cpi5fzZONSWggjjBZuj1jGbZ8s8etuXd0p3bCdAy9vIsRhozHsFNlzEnlXXYQKi+LqiekMTYoyOkRhETJK2PcCaXSuUWR0rsVkD6jgycalRNFIFK7WhT82LaVhwBzA+CtN5oxJZM6YhMPhYPOfX+LE6mrG6k8IyTjOjrZL+btjCBcMS+TiUR6u1ytX1qAga/gKq5MkalLJdSW0RJ47p8YeGUpyXQlmSKId7HY7M+5sH4x06Bjb/7yS2JWHmKwKiZgKjx+YS1xsPNdMSifJ3fV65cpqHf0o7HRevKnjY750KYZ0WQjhKUmiZpWVRRjnjsgOw9xzahKHZTDn4dsA2PnyO1R+cojJzk20xJzgvdOTKY3IZW7eQCYN7WG9XrmyWkc/CzuBNvdaBCdJomZlokn5nphw3eVwHTTXNrDh0b+R9v4JkvVhKouqeWzEfAanDmThhMFEh5/3EZQrq3l1rnVCvws7gTz3WgQPSaJmFgCT8sNjopj1gGsT8ZJN2zmw/FMmH9xCg72CV/dlUTlwGgsnDmbEwBjXA+TKak7n1zp/+tN+F3b6XU70dr+59MMLD0gSNbsAmpSfdcEksi5wDUba9Kd/Ermphgz9Hns2VPD6pPnkDRvC3LxUQi1cAw9IXTWxP/QQnN8k70Fhx+Nyorf7zaUfXnhIprgEqIoKKCx0/T8/37w56HRpGduffJXQplCaqOFkXhwNo2dzTWYUg04fl1qBGWzZ4lrxo7rTwhqxsXDvvfDLX/pk4/AeeXsRfxNtCuArgTTFJT8/P7ewsLCouLg47MMPPxxw6623nvH0uaqqqmzTp0/P7fj91KlToYsXLz7z7LPPHj3/vjLFJYgsXw7f+Y7rugauwvVzz5mzYJ2Umc7cR74PwI5X3iL0o8PY9qxhY+QpDk2fzfQExfRkLev1Gqm7JvZbbnH983cTqLf7zaUf3lLO3wqtL0nU4XBgt/877SUkJDiLior2dvw+ZsyY0dddd11VX+KRZf8CTEUFfPe7/06g4Lr+dbXvs9k2GJ/41auY/Yc7mP7wfxATncKkD/ZQ//jfeP43j/H8hsNUN7b2/iSBzoiT1tMGpD7aOLxH3u43l354S4mKisoHuP/++9O3bt06IDc3N+/BBx8c6HA4uOWWWzLGjh07etSoUXkPP/xwMrjWzp02bdqoBQsWZOfk5Izp7nk///zz8NOnT4f2dT9RqYkGmJISCOli5T2b7dyCtZm7gCKiopjz4B0AHN68DV5YT+g/3uHjf56hdMZULp05jbHpcQZHaYBOJ023tHDop8uIvWWJf/KXmQa5eXvkusVHwgcrb2yF1tnzzz+fuHDhwjM2W9/qloYmUaXU5cDvgBDgGa31r42MJxBkZUFb25dvdzr/XbA2w1RMdwdCZk+bTPa0yTgcDjb+6TkGrNvPqU+K2ZFYT/Q113PFuPTgWK/3vJOmgLT/WUruQ3P4zbMp/ikAmWmQm7eTupkKCRZSu6Esqa2q2WtboYUkhDfHzEj321Zona1YsSLxueeec3sLtA6GJVGlVAjwR2AucAzYopRaqbXe2/Mj+y6YRq6npMCzz365T7RzwdroLiBPasF2u52Z33dtIn766FEcj79E2LJVfEA9ZZOGc/niyxmSGMDr9XZx0loJJbWpxLCNCXyiL19Wbyd1MxUSLMLThOcLnmyF1mHjxo2RbW1taubMmQ19Pa6RfaIFwAGt9SGtdQvwT2CRtw9i5HZiRlmyBMrK4N13Xf+OHTs3SRnZBdS5QlVd7frZVX9tT5KGDOGyx37IpX+8jbSLh5FdeIwD//McL93zK9bsLKHNaZ0R527r4qSF0koJWf8qAFleMH5Zhce8sRVah7/97W+Jixcv9miUr5HNuelA52HEx4Bp599JKXUzcDPA0KFD+3QAMzRbGiUlBebN6/5vRnUBebsWnP+1a+Fr0Fhby/rfPo396dW8r1s4lRvPld9d0u16vZZrnWg/aXrpUmoaQwmlle+yjEpSiAyEMTD++rJa7sSL7nhrKzSAlStXJr755pv7PYnDsHmiSqnrgMu01je1/34DUKC1vrO7x/R1nmh309vWrsUU24kZzYjriT+m5B3a+An7/7GFMGckNSF1JH97EV+ZPOpf02TMPKiqVxUVvPu/JXzvoSyqw1L8Oj3Tp/r7ZXXnw2zpE98/gTRP1CjdzRP1OIkqpW7UWv/F04CUUtOBB7TWl7X/fh+A1vpX3T2mr0k0COZQW1LHtczXc/S108knj/0ex+EQtFZUDQ5hyrdvZExuhOU/EwFXoerPl9Wd5BjkFwNJov3XXRLtT5/og/14LMAWYKRSKlspFQZ8A1jZz+c8R0/T24RxlixxXbvWrnX99FVlQNlsXPyju5n9xzuZeM8sBlSe5fCvn2X5tx5n6ujN/7qfFfsUjZie6RF35rV2lAgef7zvX1Z3O9k7+hE6s+KJF6bTY5+oUmpXd38CUvtzYK21Qyl1B/Aurikuz2qt9/TnObsiI9fNyd8DIZNH5HHFE3lUlGv+v2/+Hw/O3EjYxdvYc7aVn711E1lZ0f4LJli4U0M8/z6PPw6TJrn/ZXW3k72/o+kCrurvMafT6VQ2my0AR+91z+l0KqDLEb49NucqpU4BlwHnL4OkgE+11oO9FaQ7ZO1c4Q0d1+3U2DL++7JnyI5OxkEbcRePYNrXrzQ6vMDgTvOpN5pY+/IcnvYjBEBfqhebc1cOGjQoLyUlpTpYEqnT6VQVFRVxJ0+e3DthwoSLBOP9AAAb2ElEQVSF5/+9t9G5q4ABWusd5/9BKfWRl2IUwq/+3TqRTlbWz0lJ1nz+9oscXbWbDz4poTG8la/c8zUSMtKMDtW63KkhemOodl+GmnvSLNXdqOGkJHPv7OAjDofjppMnTz5z8uTJsQTPsrFOYLfD4bipqz/KLi5CdNJWX8XqRx4n8lQKStkIzYnnK7d9/ZxFq4Ub/FUT7Xw8XzS3djVqGCA62rUMmEVqpd6qiYovkyQqRFe05otPV7Pvn1uJcabgCGkj74ZLyJg6zujIrMOd5lN/DdX2VFeJvjOLjPCVJOo7kkSF6EVz/Vne+/PjhB2KI0xFogfamXHPNwmPCeBlBr3FnRqi2QftdCR6mw3q68/9m0UmnksS9R1JokL0wZ4t71P0j7UktGSgbZB6aS5jr5ltdFjC1zp2uV+0CJqa/n271ESDniTRbpi9cCyMVVNzlo9f+D1qdyjRxOOM1ky+/WriM/06YF34m9mbn7shSdR3JIl2IQBGtAs/0VpTuHUDB159lYSabOy2UCJHJzLllmtlMFKgsmAJW5Ko70gSPY9VVgdz93tswe+7ZZ06XcXmlf9L2/ZG4p2DcIZqRn3jQoZcMN7o0ESQkyTqO8Eyz8dtVlgdzN0do2RnKf9KTUpg4Y0/Yf4TP8f+zXGcTvqCPX9dw0d3PMXHDy6jpa6bEZ5CCMuSmuh5zF4TdTc+s7+OYHHgxGl2vfc8bbtOktSUCSEhDJw5nLHXzTU6NP+QphBTkJqo70inzXn8tdemp9cWdxd58fa+ncIzI9KSGPGde2hocfDBho00b3uRxo+rqPz4IM4oRf6t80kYlm50mL4hgwtEEJCaaDd8WYDuz7VFaqL9Y3TFSGvNjpJyDm94ibaifSRUjcJuCyNyZByTb/tq4AxGkg+gqUhN1HekT7Qbvtpqyt2dm3qKy53t3WQbuC8zQx+xUor87FSu+dadzPrJI5xeNIGa3D2cPPQF6+9+hg/ufoqjG3b6PzBvs8Lggs7c2bJNiC4EfU3U3zWTrpbi9GTRE0NG5xpdjesHM1eM2pyadUXHqdi+krCju4goG0k40diS7RT84OtExA0wNkBPmPkNP18QNDtLTdR3gjqJGvHdsdK15Rx9eLPMmGu9VXjxtaNnGlizcRvpZW/QeCiE2LoslM3GwBnDGPN1iw1GssLCBJb9QvaNJFHfCdokauR3xwrXlnP04c0ya6HeatfKZkcb7+w8RnPReySd3U5zcSaRzniIVEy85SoSRmQYHaJ7zFii6swqpat+kiTqO4aMYlBKXQc8AIwGCrTWfl8Q18jRq55sa2goN9+s7rZenDPH+Nfor1HX3hJuD2HR5EyY/D32HK9m3dZd5J15m9qyBrY90YpdhRExPI4pd5h8ZaSUFPO+yeD6Ara0nHtba6vrdiHcYNS3bzdwDfC/Bh3f8O+O2a8t53DzzTL7tBojCi/eqIiNGRzHmIUzqWm6gDe2HyGs5GOGt27nxO5k1t39DNhh+DUFDL1okjdDDw6elq7MXsMWfmNoc65S6iPgR+7WRH3VJ2qZZlUjufFmWa3J1Nd81bSttWbjodNs3V1EQc0ammuraPoimwiiUQl2Cn5wLREJcf0/UDDpS1I0a59FD6Q513dMn0SVUjcDNwMMHTp0cmlpqVdjkAJlH7jxZknBxMVfBYqT1U2s2H6UhPJNTNLbKP4iiujTQ7DZbCRNy2LcNy/z3sGEZUuKkkR9x2fNuUqptcCgLv50v9b6DXefR2v9NPA0uGqiXgrvXyzVrGo0N94sy/X3+oi/mrYHxUXwn5eOpLVtOGv3zuFg5GEubvoAW2s5R7ZWc3pTKTpCM37pZSTlZnnvwMHK7H0Wwu98lkS11nN89dzC3KRg4v8+99AQG1eMS+OKcWkcKB/Hyh3HyZixnamhO9m+X1P4x1XYCSMyK5bJt1+DPSKs9ycVX2b0YAphOiYe1ieEdRk5GnjEwBjumZdDQ8twVu6YwQmOc9mITxjAQXZsSmLdvX+BEMheMIms2YEzjcMvrDbMW/icIX2iSqnFwJNACnAW2KG17rXzxp9r5wrhDWboc9das/3IWdbuPcmolr3MsRey42QT9TsziNRRqHg7U+66lqhkGYzkNjOc2D6QPlHfCdrFFoQIRqfrmllRWEbt2QoWhWwkru0Y63cOIOZUKjZlI2lqJuNuuNzoMIWXSRL1HWnO7YP+FD4tVnAVASppQDg3zRxGmzObj/dlsengafInH2CGbTtFNQ2Ubqrl9GdHIBzG/scckscOMzpkIUxNaqJu6s/UMAtOKxNB5OiZBl4vLEM31fDViM9Ichxl7X47YcUDCSOMiCExTP7+tTIYycKkJuo7kkTd0J+pYRadViaCUFNrG2/vPkHRiVqmRx3lwtZNlLU2UbghlriaBFQIZF2ZT9a8AqNDFX0kSdR3pDnXDf2ZGibTyoRVRISGsDg/A/Jhd9lgntidSUxoE9ddvoW4hhLWnQnni9WfUvrmDlScncl3XE30oESjwxbCUJJE3ZCVBc3N597m7tQwmVYm/Mlbfe9j0+MYmx5HdWMrrxfGcbytkdnZ5Vw8cB1VzjN8XBjJpodexK5CSMzPIO/blxESEuKtlyGEZUgSdcPateB0/vv30FD3p4bJtDL/C9ZBXL7oe4+LDOU7X8lyrdd7MIVfVw1kUKSTr128jaiz+yhsi+bwh7WcvvsYKkwx+oZLGDhxpHdekBAWIH2ivfBWn2awXtj9LVgHcfmz7/1kdROvFR6jtsnBwrRqcs+8T2NbK2uPRBOxI4JwHUZExgDyv38NYZER3j248Ij0ifqO1ER74a0+TX8theerZO3PQoCnxzLzfqa+5s++90FxEdx2yQha25ys2XuKFY2LGZloZ9HEbTiS9rCtPoEje6tp+PFzqBAbQ+eNZ9iVF3g3CCFMwmZ0AGZnpT7N5ctdtZG5c10/ly839/N6+1gdiaSzjkQS6Iz4nIaG2LhyXBo/vXI0+cMGcevHk8h84gZ++KtZFG1ykjjyFI7ZDore28THdz7NuvufpfZ4pe8CEsIA0pzrBits7+Wr5jx/NhP291jBPp3IyM9px3vf5HAQnXcce0I9oaejKf7zZyTU76QyZiAfb24h/mAUdmUnflwaY268HLtdGsP8QZpzfUc+wW6wwvZevmrO82czYX+PFeyDuIz8nP773Nmp2zkU0CSMrOKB4rGMHp7PdYNbuDZ/FXrMKbYOyKL09WKqfnAcW5gi5/qLSJ2c479gzUIGSgQEqYkGCKmJnvs8cm3yr57OnS3StV5vRV0zs0bEU9CyGXV8O/UDUllzGKI3acJ0OBFp0eTfdQ1h0UEwGMnPI+CkJuo7kkQDiK+a8/zZTNjXY0nCNI/ezl2bU/PxvnI2HzpDRmIUi7McDNi3AlrqKE7JY/cbpSRWxGILsZFx6WiGL7rQuBfjSwb0O0gS9R1JogHGzKNz3X0Od+8XrNNZzMzdc3fkdAOv7yijtc3JlWNSGF27CY5uonVAKu/Xx+J8q5JoRzS2AXYm/Od8YocO9NdL8L0tW1wj56qr/31bbKxrQvpU3+zvKknUdySJ9oHUejzn7YTXY2EeOVFW0dTaxurPT1B8spa8wbFcPsRB+N5Xoamak0Mms27NQRL3hWNXocSNTmXs9660/mAkqYkGFEmibpJaj+d8cc3orjC//d7lDP+lnCgr2l1WzTu7TxIRauPqCYPIOL0RStfjjB7IpuhMypcXE18fiwpVjPrGDNIK8owO2XN+HkotSdR3DEmiSqmHgQVAC3AQuFFrfba3x1lxFxfhm9arrs7JkIgKSlUmSk6UpbnW6y3j+NlGpg9P4qJBDmy7X4aGM9QMLWDNrmPErG8hXEcQnhrFxDuvJiI22uiw+86PTVuSRH3HqMUW1gBjtdbjgX3AfQbF4ZZgnsTvDb5YCKBjOktkpCshR0bC/91fgpITZXkd6/X+5IpcQkNs/ObTapbphVRNv4/YkAiujS9j3tc0yT/I4RjH2Xj/3/jkrv9j/2ufGB1636SkuEqRUsCzNMObc5VSi4Gvaq2/2dt9pSbaO7P22/qq9eqc14uFTpTokxPVjby2vYy6ZgeXjRnEhPhm1OcvQ305zUOn8355Daw4TnRrNCHRdsbdcgVx2WlGh20aUhP1HTMk0TeBF7XWf+/m7zcDNwMMHTp0cmlpqT/D+xcrrFpk9n5bvyR4L58osxZKLMEHb15rm5P39pxi57GzjEgZwILxg4gs+xQOfgCRCZQNv4h1r39Myt5wQlUYsaOSGXvrfOsPRuonSaK+47MkqpRaCwzq4k/3a63faL/P/cAU4BrtRiAyOrd77tSWzRy/V3nphZq9UGJqfnjz9p+q5c2dx1FKsWjiYIZFNcPnL0HNcdqGXsCnhHP6+ULi6+KwhdkYfu0FpH9lrFdjsApJor5jWE1UKfUd4FZgtta6wZ3HGJ1Ezay3wTuSEPrGSk34puPnN6++2cEbO45TcrqeyZkJzM5JwX5sExxYA+ExVI26jDXr1hO3rpVIHUFYciQTv7+QiPhYr8diVpJEfceo0bmXA48BF2utK9x9nCTR7vV03QJJCH1lwHz4wGHQm6e1ZvuRKtZ+UU58ZCiLJ6UzMKQBPn8Fqo+iMwr4PC6Nouc/IOVELCE2O4MuHMGo6y7xWUxmIUnUd4zqKPgDEA6sUUoBbNJa32pQLAGhp8XXt2zx3yLygcJKW+CZjkFvnlKKyZmJTM5MpLKumRXbyzhd38KlOddSMDUBdWwL4794i/EzomgcfQnv7Suk9rVCTnyyn5DoEPKWziNxZIZPYxSBx/CBRX0hNdHeddUdKE2TnrHCYDLTMsmb1+bUfFRczmeHzzAkMYpFEwcTo+th9ytQVQrpkyhNG8enr7xByu5IwggjZkQi425bGFCDkaQm6juSRIOESa5plhM0g7F8wWRvXunpet7YcZzWNidXjU8jNzUGjm+HotVgD8eRdzXrzhyh6vlCEmtjsYWGMGxRARkXjzc69H6TJOo7kkSDiMmuaUFH3n9zaGpt461dJyg+VcuYwbFcMTaNsLZ62P0qnD4IaROozPoKa95/ncQPHUQ6IwlNjmDi7QuJTI4zOnyPSBL1HUmi3ZALnvAmGR1tTp8fq+adPSeIDA1h8aQM0uMj4fgOKFoFNjs672q2O2o58Nf3GVgWT0iInZSCbEZfP9vo0PtEkqjvSBLtQqBd8KRAYCzpkza/6oZWVhQe43h1EzNGJDNzRDI2RwPsWQEVxZA6hvqRc3i38D1CV5wktjkGW2QIY5bOJjEn0+jweyVJ1HckiZ4n0C54gVYgsCKZLmMdTqfm04OnWbe/goGxEVw7KZ34qDA4uRv2vuG605irORgewcaXX2fQrgjCbOFEZ8Yz4fZF2MNDjX0B3ZAk6juSRM8TSBe8QCsQWJWcB2vqWK+3vmO93iHx0NLgSqbleyElh9bc+XxYupG6v39OYnU8IaEhZM6fyNBZk40O/xySRH0ncMZwe0kgzQ/s2H1G5ocaq6c5vMK80uIiuf3SEbQ4nLy39ySrPz/B8IEDWDD+a0RODIHyIkI3P8U8pwNuX8Cp2EGsWfsK9W+sp+T1QkITwhl/5wKik+ONfinCh6Qm2oVAmQ4iNSBzkb5p69t3qpZVu05gU7BoYjrZydHQ2uQaiHRiJyQNxzlmMduqijnw9w9IPRJPaIidpMlDybthnmFxS03UdySJdiNQLniBUiAQwkzqmh2sbF+vd0pmArNyB2IPsUHlAddUmbZmyL2KmuQRvLP9DSJfrSC2JZaQiBByv3MpyWOy/RqvJFHfkSQaBIwoEARKISSYyTnsndaabaVVvF/Uab3emAhwtEDxW1C2HRKyYNxXKW44yeYVb5K2M5JwFU7U0DjG376QsIhwn8cpSdR3JIkKr5MRwdYn57DvKmqbeb3QtV7vrNyBTM1KQCkFZw65aqctDZBzBc1p43n/wPs0/GMPKVXx2OwhDLliHFnzCnwWmyRR35EkalFmrSVIP6z1yTnsnzan5sOicraUdFqvNyIU2lph3ztw9DOIHwrjvsrxtkbWfPAKqR9oopzRhMaFMf72+UQPSvRqTJJEfUdG51qQmWsJRo8INqpwYdZCjSeMPodWF2JTzMlLZU5eKiWV9Ty7vgSH08n88YPJGb0ARi+As0dg67MMbq7lOxPm0XZVAZtPbObQP9fR+NBLhNpCSZg4mNHfnkdISIjRL0n0QGqiFmP2WoKR8RlVuDBzocYTZv+MWVFTaxurdp1g36laxqbHcfmYQYTZbdDmcG0efmQjxKTBuK9xNiSE1dtXMGDFaWKaY7jw0e8SGtq/RRykJuo7kkQtxgqLQRgxItioC3+gJhwrj+o2e6vArmNneWf3SaLCOq3XC1BdBp+/DI1VMGI2OvNCDteUMCx+WL+PKUnUd6Q512KssBjEkiUwZ45/L2RGNUEGatOnEefQG6zQKjA+I57xGfFUN7TyWuExTrav13vhiMHYLrwbnG1w8APU+w8wLCYNpt0KShkdtuiGITVRpdQvgEWAEygH/kNrfby3x0lN1MXKtQRfkZqosOq5cDo1Gw5Wsv5AJakxEVzTsV4vQONZiOz/ikdSE/Udm0HHfVhrPV5rPRFYBfzMoDgsackS14Vh7VrXT7Mk0IoKV3NzRYX/j92xtF5kpKt5OzLSx0vrtb/YFCr8e1zRrY5Wgc46WgXMzGZTzByZwn1XjObysYN4YfMRfvtOEbuOeSeBCt8yvE9UKXUfMFRr/Z+93VdqouZllmY0v/SHdfFiK+YssVzTZ6Cxak20Kx3r9R4or+Ou2SNd8037QWqivmNYElVKPQR8G6gGLtVad1l/UUrdDNwMMHTo0MmlpaX+C1K4JZAuXr0KqhdrPdLV0TVJor7js+ZcpdRapdTuLv4tAtBa36+1HgK8ANzR3fNorZ/WWk/RWk9JkYuUKVm1Gc0jQfVircesXR0icPlsdK7Weo6bd/0H8Bbwc1/FInzLCiOGvSaoXqw1paRIo4DwH0MGFimlRnb6dSFQZEQcwjv8PqjHSEH1YoUQvTFqisurQA6uKS6lwK1a67LeHicDi8zN7JPcvSqoXqywOukT9R1DFlvQWl9rxHGFbwVVM1pQvVghRHeMmicqhBBCWJ4kUSEsxshFLYQQ55IkKoSFLF/umqY6d67r5/LlRkckRHCTJCpMyWq1LX/EW1HhWkigsdG1i09jo+t3q7xHQgQiSaLCdKxW2/JXvLLOgxDmY/jauX0hU1wCn9VW1fNnvFZ7b4R5yBQX35GaqDAVq9W2/BmvrPMghPnIptzCVKy2qp6/47XqZtlCBCqpiQrT+elPISLCGrUtI2qHKSkwdap53xMhgonURIVpdN6mUym491645RbzJwupHQoRvGRgkTAFGTQjhO/IwCLfkeZcYQpWG1AkhBAgSVSYhNUGFAkhBEgSFSYh0zeEEFYkA4uEacgAHSGE1UgSFaYi23QKIazE0OZcpdSPlFJaKZVsZBxCCCGEJwxLokqpIcBc4IhRMQj/s9ruLEII0RMja6KPAz8GrDNRVfSL1XZnEUKI3hiSRJVSC4EyrfVON+57s1Jqq1Jqa4VUXyxL9sIUQgQinw0sUkqtBQZ18af7gZ8C89x5Hq3108DT4FqxyGsBCr/qWEyh84pEHYspyEAiIYRV+SyJaq3ndHW7UmockA3sVEoBZADblVIFWuuTvopHGEsWUxBCBCK/N+dqrT/XWg/UWmdprbOAY8AkSaCBTRZTEEIEIpknKvxGFlMQQgQaw5Noe21UBAlZTEEIEUhk7VwhhBDCQ5JEhRBCCA9JEhVCCCE8JElUCCGE8JAkUSGEEMJDSmvrLAKklKoASo2OA0gGKo0Owg0Sp/dZJVaJ0/usEmtXcWZqrWVcvA9YKomahVJqq9Z6itFx9Ebi9D6rxCpxep9VYrVKnIFCmnOFEEIID0kSFUIIITwkSdQzTxsdgJskTu+zSqwSp/dZJVarxBkQpE9UCCGE8JDURIUQQggPSRLtA6XUs0qpcqXUbqNj6YlSaohS6kOl1BdKqT1KqbuMjqkrSqkIpdRnSqmd7XE+aHRMPVFKhSilCpVSq4yOpSdKqRKl1OdKqR1Kqa1Gx9MdpVS8UuoVpVRR+2d1utExnU8pldP+Pnb8q1FK3W10XN1RSv2g/bu0Wym1XCkVYXRMgU6ac/tAKXURUAf8VWs91uh4uqOUSgPStNbblVIxwDbgaq31XoNDO4dy7coerbWuU0qFAuuBu7TWmwwOrUtKqXuAKUCs1nq+0fF0RylVAkzRWpt6TqNS6nlgndb6GaVUGBCltT5rdFzdUUqFAGXANK21Gearn0MplY7rO5SntW5USr0ErNZaP2dsZIFNaqJ9oLX+BDhjdBy90Vqf0Fpvb/9/LfAFkG5sVF+mXerafw1t/2fKUp1SKgO4CnjG6FgCgVIqFrgIWAagtW4xcwJtNxs4aMYE2okdiFRK2YEo4LjB8QQ8SaIBTimVBeQDm42NpGvtTaQ7gHJgjdbalHECTwA/BpxGB+IGDbynlNqmlLrZ6GC6MQyoAP7S3kT+jFIq2uigevENYLnRQXRHa10GPAIcAU4A1Vrr94yNKvBJEg1gSqkBwKvA3VrrGqPj6YrWuk1rPRHIAAqUUqZrJldKzQfKtdbbjI7FTTO01pOAK4Db27shzMYOTAL+rLXOB+qBnxgbUvfam5sXAi8bHUt3lFIJwCIgGxgMRCulvmVsVIFPkmiAau9jfBV4QWv9mtHx9Ka9Ke8j4HKDQ+nKDGBhe1/jP4FZSqm/GxtS97TWx9t/lgMrgAJjI+rSMeBYp5aHV3AlVbO6AtiutT5ldCA9mAMc1lpXaK1bgdeArxgcU8CTJBqA2gfsLAO+0Fo/ZnQ83VFKpSil4tv/H4nrIlBkbFRfprW+T2udobXOwtWk94HW2pQlfKVUdPtgMtqbR+cBphtNrrU+CRxVSuW03zQbMNXAt/MswcRNue2OABcopaLarwGzcY2HED4kSbQPlFLLgY1AjlLqmFJqqdExdWMGcAOuGlPH0PwrjQ6qC2nAh0qpXcAWXH2ipp4+YgGpwHql1E7gM+AtrfU7BsfUnTuBF9rP/0TglwbH0yWlVBQwF1fNzrTaa/WvANuBz3Fd32X1Ih+TKS5CCCGEh6QmKoQQQnhIkqgQQgjhIUmiQgghhIckiQohhBAekiQqhBBCeEiSqBD9pJS6XClVrJQ6oJQy7ao7QgjvkykuQvRD+84e+3DNIzyGa77rErPtmCOE8A2piQrRPwXAAa31Ia11C65lARcZHJMQwk8kiQrRP+nA0U6/H8OE284JIXxDkqgQ/aO6uE36SIQIEpJEheifY8CQTr9nIBshCxE0JIkK0T9bgJFKqez2PSe/Aaw0OCYhhJ/YjQ5ACCvTWjuUUncA7wIhwLNa6z0GhyWE8BOZ4iKEEEJ4SJpzhRBCCA9JEhVCCCE8JElUCCGE8JAkUSGEEMJDkkSFEEIID0kSFUIIITwkSVQIIYTwkCRRIYQQwkP/P6aO9pttYF7YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.axes()\n",
    "\n",
    "df_X.query('label == -1').plot.scatter(x=0, y=1, ax=ax, color='blue')\n",
    "df_X.query('label == 1').plot.scatter(x=0, y=1, ax=ax, color='red')\n",
    "\n",
    "_xs = np.array([np.min(Xs[:,1]), np.max(Xs[:,1])])\n",
    "for k, theta in enumerate(all_thetas):\n",
    "    _ys = (theta[0] + theta[1] * _xs) / (- theta[2])\n",
    "    plt.plot(_xs, _ys, label='iter {0}'.format(k + 1), lw=0.5)\n",
    "plt.legend(bbox_to_anchor=(1.04,1), loc=\"upper left\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since all the lines are so close to each other, Newton's method converges very quickly, at least for this dataset."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
